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Abstract 

We present a general methodology for constructing lattice Boltzmann models of hydrody- 
namics with certain desired features of statistical physics and kinetic theory. We show how a 
methodology of linear programming theory, known as Fourier-Motzkin elimination, provides an 
important tool for visualizing the state space of lattice Boltzmann algorithms that conserve a 
given set of moments of the distribution function. We show how such models can be endowed 
with a Lyapunov functional, analogous to Boltzmann's H , resulting in unconditional numeri- 
cal stability. Using the Chapman-Enskog analysis and numerical simulation, we demonstrate 
that such entropically stabilized lattice Boltzmann algorithms, while fully explicit and perfectly 
conservative, may achieve remarkably low values for transport coefficients, such as viscosity. 
Indeed, the lowest such attainable values are limited only by considerations of accuracy, rather 
than stability. The method thus holds promise for high-Reynolds number simulations of the 
Navier-Stokes equations. 



Keywords: computational fluid dynamics, thermodynamics, hydrodynamics, entropy, numerical 
stability, lattice gases, lattice Boltzmann equation, detailed balance, complex fluids 
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1 Introduction 



Though lattice models have been used to study equilibrium systems since the 1920's, their ap- 
plication to hydrodynamic systems is a much more recent phenomenon. Lattice-gas models of 
hydrodynamics were first advanced in the late 1980's, and the related lattice Boltzmann method 
was developed in the early 1990's. This paper describes a particularly interesting and useful cate- 
gory of lattice Boltzmann models, which we call Entropic Lattice Boltzmann models, and provides 
a number of mathematical tools for constructing and studying them. This introductory section de- 
scribes and contrasts lattice methods for hydrodynamics, traces their development, and endeavors 
to place the current work in historical context. 

1.1 Lattice Gases 

The first isotropic lattice-gas models of hydrodynamics were introduced in the late 1980's 
Such models consist of discrete particles moving and colliding on a lattice, conserving mass and 
momentum as they do. If the lattice has sufficient symmetry, it can be shown that the density 
and hydrodynamic velocity of the particles satisfy the Navier-Stokes equations in the appropriate 
scaling limit. 

This method exploits an interesting hypothesis of kinetic theory: The Navier-Stokes equations 
are the dynamic renormalization group fixed point for the hydrodynamic behavior of a system of 
particles whose collisions conserve mass and momentum. We refer to this assertion as a hypothesis, 
because it is notoriously difficult to demonstrate in any rigorous fashion. Nevertheless, it is very 
compelling because it explains why a wide range of real fluids with dramatically different molecular 
properties - such as air, water, honey and oil - can all be described by the Navier-Stokes equations. 
A lattice gas can then be understood as a "minimalist" construction of such a set of interacting 
particles. Viewed in this way, it is perhaps less surprising that they satisfy the Navier-Stokes 
equations in the scaling limit. 

For a time, lattice-gas models were actively investigated as an alternative methodology for com- 
putational fluid dynamics (CFD) 0j. Unlike all prior CFD methodologies, they do not begin with 
the Navier-Stokes equations; rather, these equations are an emergent property of the particulate 
model. Consequently, the use of such models to simulate the Navier-Stokes equations tends to 
parallel theoretical and experimental studies of natural fluids. One derives a Boltzmann equation 
for the lattice-gas particles, and applies the Chapman-Enskog analysis to determine the form of the 
hydrodynamic equations and the transport coefficients. In fact, numerical experimentalists often 
circumvent this theory by simply measuring lattice-gas transport coefficients in a controlled setting 
in advance of using them to study a particular flow problem, just as do laboratory experimentalists. 

One often overlooked advantage of lattice-gas models is their unconditional stability. By insist- 
ing that lattice-gas collisions obey a detailed-balance condition 0, we are ensured of the validity 
of Boltzmann's H theorem, the fluctuation-dissipation theorem, Onsager reciprocity, and a host 
of other critically important properties with macroscopic consequences. By contrast, when the 
microscopic origins of the Navier-Stokes equations are cavalierly ignored, and they are "chopped 
up" into finite-difference schemes, these important properties can be lost. The discretized evolution 
equations need no longer possess a Lyapunov functional, and the notion of underlying fluctuations 
may lose meaning altogether. As the first practioneers of finite-difference methods found in the 
1940's and 1950's, the result can be the development and growth of high-wavenumber numerical 

x This is certainly possible for an ideal single-phase viscous fluid. For complex fluids with finite-range interaction 
potentials, it is an outstanding problem. 
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instabilities, and indeed these have plagued essentially all CFD methodologies in all of the decades 
since. Such instabilities are entirely unphysical because they indicate the absence of a Lyapunov 
functional, analogous to Boltzmann's H. Numerical analyists have responded to this problem with 
methods for mitigating these anomalies - such as upwind differencing - but from a physical point 
of view it would have been much better if the original discretization process had retained more of 
the underlying kinetics, so that these problems had not occurred in the first place. Lattice gases 
represent an important first step in this direction. As was shown shortly after their first appli- 
cations to hydrodynamics, lattice-gas models for single-phase ideal fluids can be constructed with 
an H theorem that rigorously precludes any kind of numerical instability. More glibly stated, 
lattice gases avoid numerical instabilities in precisely the same way that Nature herself does so. 

1.2 Lattice Boltzmann Methods 

In spite of these appealing features, the presence of intrinsic kinetic fluctuations makes lattice- gas 
models less than ideal as a CFD methodology. Accurate values for the velocity field at selected 
locations, or even for bulk coefficients such as drag and lift, have an intrinsic statistical error that 
can be reduced only by extensive averaging. On the other hand, the presence of such fluctuations 
in a simple hydrodynamic model make lattice gases an ideal tool for those studying the statistical 
physics of fluids, molecular hydrodynamics, and complex-fluid hydrodynamics; not surprisingly, 
these remain the method's principal application areas Consequently, many CFD researchers who 
appreciated the emergent nature of lattice-gas hydrodynamics but wanted to eliminate (or at least 
control) the level of fluctuations, turned their attention to the direct simulation of the Boltzmann 
equation for lattice gases. Such simulations are called lattice Boltzmann models (71 l8l l9l FTTH \T\] . 
These models evolve the real-valued single-particle distribution function, rather than the discrete 
particles themselves, and in this manner eliminate the kinetic fluctuations. 

Early attempts along these lines restricted attention to Boltzmann equations for actual lattice 
gases. It was soon realized, however, that these quickly become unwieldy as the number of possible 
particle velocities increases. In order for lattice Boltzmann models to become practical tools, it 
was necessary to develop simplified collision operators that did not necessarily correspond to an 
underlying lattice-gas model. The most successful collision operators of this type are those of the 
form due to Bhatnager, Gross and Krook [T2|, and these have given rise to the so-called lattice 
BGK models H3J. 

1.3 Lattice BGK Models 

Lattice BGK collision operators allow the user to specify the form of the equilibrium distribution 
function to which the fluid should relax. For lattice gases obeying a detailed balance H condition, 
this is known to be a Fermi-Dirac distribution. Having abandoned lattice-gas collision operators, 
however, it seemed unnecessary to continue to use lattice-gas equilibria, and practitioners exploited 
the freedom of choosing lattice BGK equilibria to achieve certain desiderata, such as Galilean 
invariance, and the correct form of the compressible Navier-Stokes equations. 

The alert reader will have noticed, however, that the evolution to lattice-BGK methods has 
jettisoned the last vestiges of kinetic-level physics that might have been left over from the original 
lattice-gas models. The move to a Boltzmann description of the lattice gas eliminated kinetic 
fluctuations, but at least it retained an i?-theorem; more precisely, its global equilibrium still 
extremized a Lyapunov functional of the dynamics. The move to lattice-BGK operators and the 

2 Actually, a weaker condition called semi- detailed balance suffices for this purpose. 
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arbitrary legislation of the equilibrium distribution function, however, completely abandoned even 
the concept of detailed balance and the ff-theorem. Without a Lyapunov functional, lattice-BGK 
methods became susceptible to a wide variety of numerical instabilities which are ill-understood 
and remain the principal obstacle to the wider application of the technique to the present day. 

1.4 Motivation for Entropic Lattice Boltzmann Methods 

In this paper, we argue that the most natural way to eliminate numerical instabilities from lattice 
Boltzmann models is to return to the method its kinetic underpinnings, including the notion of a 
Lyapunov functional. We present a general program for the construction of "entropically stabilized" 
lattice Boltzmann models, and illustrate its application to several example problems. 

The crux of the idea is to encourage the model builder to specify an appropriate Lyapunov func- 
tional for the model, rather than try to blindly dictate an equilibrium state. Of course, specifying 
a Lyapunov functional determines the equilibrium distribution, but it also governs the approach to 
this equilibrium. It can therefore be used to control the stability properties of the model. It should 
be emphasized that the presence of a Lyapunov functional guarantees the nonlinear stability of the 
model, which is a much stronger condition than linear stability. 

As we shall show, collision operators that are constructed in this way are similar in form to the 
lattice BGK collision operators, except that their relaxation parameter may depend on the current 
state. As a result, the transport coefficients may have a certain minimum value in models of this 
type but, as we shall see, this value may actually be zero. In fact, the ultimate limitation to the 
application of these algorithms for very small transport coefficients will come from considerations 
of accuracy, rather than stability. 

1.5 Structure of Paper 

The plan of this paper is as follows: In Sectional we present the general program of construction of 
entropic lattice Boltzmann models. This includes the introduction of a conservation representation 
of the lattice Boltzmann distribution function, by means of which the collision process is most 
naturally described. This representation also provides an interesting geometric interpretation, as it 
allows us to describe a collision process as a mathematical map from a certain polytope into (and 
perhaps onto) itself. We then show a way of constructing Lyapunov functionals on this polytope, 
contrast these with Boltzmann's H, and use them to construct collision operators for absolutely 
stable models. 

Section |H1 carries out this program for a very simple lattice Boltzmann model of diffusion in one 
dimension. The model is useful from a pedagogical point of view, since an entropically stabilized 
collision operator can be written for it in closed form. Moreover, it is possible to determine the 
transport coefficient of this model analytically. The result is a fully explicit, perfectly conservative, 
absolutely stable method of integrating the diffusion equation. We present numerical simulations 
for various values of the diffusivity in order to probe the limitations of the technique. In the course 
of presenting this model, we show that a linear algebraic procedure, called the Fourier-Motzkin 
algorithm, is very useful for constructing and visualizing the conservation representation. 

Section 0] applies the method to a simple five-velocity model of fluid dynamics in one dimension, 
first considered by Renda et al. in 1997 Here the geometric picture involves a four dimensional 
polytope, and the Fourier-Motzkin algorithm is shown to be very useful in describing it. Indeed, 
Appendix lAl gives the reader a "tour" of this polytope. 



5 



Finally, Section [5] applies the method to the two-dimensional FHP model and the three- 
dimensional FCHC model of fluid dynamics. For these examples, we are able to provide appropriate 
conservation representations. Though the latter example has too many degrees of freedom to allow 
geometric visualization of the collision dynamics, it is nevertheless possible to develop and present 
a Lyapunov functional and an entropic model. We use these more complicated models to demon- 
strate that entropic lattice Boltzmann models are not computationally prohibitive, even when the 
number of velocities involved is large. We describe the computational aspects of such models in 
some detail, and end with a discussion of Galilean invariance. 
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2 Entropic Lattice Boltzmann Models 



In this section, we present the general program of construction of entropic lattice Boltzmann models. 
We note that the lattice Boltzmann distribution function admits a variety of representations, and we 
show how to transform between these. In particular, we introduce a conservation representation, by 
means of which the collision process has a very natural description. We also examine the geometric 
aspects of this representation, and use these insights to show how to construct lattice Boltzmann 
collision operators possessing a Lyapunov functional. 

2.1 Representation 

Lattice Boltzmann models are constructed on a O dimensional regular lattice C, and evolve in 
discrete time intervals At. We imagine that a population of particles exists at each site x S C, and 
that these particles can have one of a finite number n of velocities, Cj/At where j £ {1, . . . , n}. The 
displacement vectors Cj are required to be combinations of lattice vectors with integer coefficients, 
so that if x 6 £, then x + Cj € C. We note that there is nothing precluding some of the Cj from 
being the same, or from being zero; indeed, this latter choice is often made to incorporate so-called 
"rest particles" into the model. 

Mathematically, the state of the system of particles is represented by n real numbers at each 
site x 6 C and time step t. These are denoted by Nj(x,t) £ !R, and can be thought of as the 
expected number of particles with mass mj and velocity Cj/At at site x and time t. As such, 
we expect them to be nonnegative, Nj(x,t) > 0. Taken together, these quantities constitute the 
single-particle distribution function of the system. This is all the information that is retained in a 
Boltzmann description of interacting particles. 

In what follows, we shall think of the set of single-particle distribution functions as a vector 
space. That is, we shall regard the n quantities Nj(x,t), for i £ {1, . . . ,n}, as the components of 
an n-vector or "ket," |N(x, t)) S 5i n . If no ambiguity is likely to result, we shall often omit the 
explicit dependence on t and write simply |N(x)). If discussion is restricted to a single lattice site, 
we may further abbreviate this as |N). 

In one discrete time step At, the state of the system is modified in a manner that is intended to 
model collisions between the particles at each site, followed by propagation of the collided particles 
to new sites. The collisions are modelled by modifications to the distribution function 



Here and henceforth, we use a prime to denote the postcollision state. It is usually the case that 
|N'(x)) depends only on |N(x)), though in more sophisticated lattice Boltzmann models it may also 
depend on the |N)'s in a local neighborhood of sites about x. Following convention, we define the 
collision operator Cj(|N(x))) as the difference between the old and new values of the single-particle 
distribution, 



N(x))^|N'>. 





(2) 



As this is the difference of two kets, it can also be thought of as a ket 




(3) 
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2.2 Conservation Laws 



The collision process specified in Eq. is usually required to conserve some number of locally 
denned quantities. Usually, these quantities are additive over the lattice sites and directions. For 
example, we may require that there be a conserved mass 

M=5>(x,t), (4) 

where we have defined the mass density 

n n 

p(x, i) = £ mjNj^, t) = J2 rnjN>(x, t). (5) 

j=0 j=0 

To streamline the notation for this, we can define a covector or "bra" by the prescription 

(p\ j = m j , (6) 

and write simply 

p=(p\N) = (p\-N'). (7) 

Thus, to the extent that we think of the single-particle distribution function as a vector, to each 
conserved quantity there corresponds a covector such that the value of the conserved quantity is 
the contraction of the two. Because collisions are required to preserve this value, the contraction 
of the covector with the collision operator must vanish, 

(p\C) = 0, (8) 

as can be seen from Eqs. (j3J) and Q. 

Lattice Boltzmann models may also conserve momentum, 

n n 

tt(x, t) = £ mj^-Njix, t),= J2 mj^Nfa *)■ (9) 

3=0 j=0 

This is still linear in the Nj, so that we can define 

H = m ixt> (10) 

and write 

7r(x,t) = (it I N) = <tt | N') . (11) 

Here we must be a bit careful to define the meaning of the quantities involved: Note that (tv\- is 
a covector in its j index, but a vector in its spatial index Thus, the contraction with |N)^ is 
over the j index, and results in the vector 7r. Once again, this covector annihilates the collision 
operator, 

<7T | C) = 0, (12) 

where the right-hand side is a null vector. 

3 The pedant will note that momentum is more properly thought of as a covector in its spatial index, but we do 
not bother to distinguish between spatial vectors and covectors in this paper. 
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The subsequent propagation process is described mathematically by the prescription 



Nj(x + cj,t + At) = 7Vj(x, t). 



(13) 



Note that this operation must be carried out simultaneously over the entire lattice. This may 
alter the values of the conserved quantities at each site, but because it is nothing more than a 
permutation of the values of the distribution function about the lattice, it clearly leaves unaltered 
the global values of the conserved quantities, 



Combining Eqs. © and (|13|) . we find the general dynamical equation for a lattice Boltzmann model, 



For compressible fluids, it is also necessary to pay attention to conservation of energy. Conser- 
vation of kinetic energy can be expressed using the bra, 



Indeed, the right-hand side could be generalized without difficulty to anything that depends only 
on j. The problem of including a potential energy function between particles at different sites, 
however, is much more difficult. If we suppose that there is a potential ^(|xj — x^l) between 
particles at Xj,Xfe £ C, then two problems arise: First, the total potential energy will depend 
on the pair distribution function - that of finding two particles a certain distance apart. This is 
outside the scope of Boltzmann methods, which retain only the single-particle distribution. We 
can avoid this problem by making the mean-field approximation, in which the probability of having 
one particle at Xj and another at x^ is simply the product of the two single-particle probabilities, 
but then the potential energy is quadratic, rather than linear, in the Nj's. Second, and perhaps 
more distressingly, the potential energy is not preserved by the propagation step 23 • These 
considerations make it very difficult to add potential interaction to lattice Boltzmann models. If 
we are willing to give up the idea of an exactly conserved energy, and instead consider isothermal 
systems, then methods of skirting these difficulties have been known and actively investigated for 
some time now [161117) . More recently, a method of incorporating interaction that maintains exact 
(kinetic plus potential) energy conservation has also been proposed ^Bj- I n any case, this paper 
will be restricted to systems with kinetic energy only. The possibility of extending our methods to 
models with nontrivial interaction potentials is discussed briefly in the Conclusions section. 

2.3 Geometric Viewpoint 

The requirements of maintaining the conserved quantities and the nonnegativity of the distribution 
function place very important constraints on the collision process. Much of this paper will be de- 
voted to describing - algebraically and geometrically - the most general set of collisional alterations 
of the Nj(x.) that meet these requirements. 

Suppose that there are n c < n independent conserved quantities. As described above, these 
define n c linearly independent covectors or bras, whose contraction with the collision operator 
vanishes. We shall sometimes adopt special names for these; for example, that corresponding to 
mass conservation shall be called (p|, that for momentum conservation shall be called {tt\, etc. 
Generically, however, let us refer to them as (A CT |, where a = 1, . . . ,n c . These covectors are not 





(14) 




(15) 
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necessarily orthogonal in any sense. Note, for example, that the covectors for mass and kinetic 
energy in Eqs. © and Q16JI are not orthogonal with respect to the Euclidean metric. In fact, there 
is no reason to insist on any kind of metric in this covector space. 

The n c covectors are not uniquely defined, since a linear combination of two conserved quantities 
is also conserved. The n c dimensional subspace defined by the covectors, on the other hand, is 
uniquely defined, and we shall call it the hydrodynamic subspace^. It is always possible to construct 
n — n c more covectors that are linearly independent of each other, and of the n c corresponding to 
the conserved quantities. For example, this may be done by the Gram-Schmidt procedure. Let us 
generically call these (ckJ, where n = 1, . . . ,n — n c ; if there are only a few of these and we want 
to avoid excessive subscripting, we shall use the successive Greek letters (a\ , ((3\ , (7) , . . . for this 
purpose. Once again, these are not uniquely defined or orthogonal in any sense, but they do span 
an n — n c dimensional subspace which we shall call the kinetic subspace. The union of all n of 
our covectors, namely the (A^j's and the (cell's, thus constitute a basis for the full n dimensional 
covector space. 

The ket |N) whose components are the single-particle distribution function is thus seen to be 
but one representation of the dynamical variable. As has been recognized for some time now |19j . 
we are free to make a change of basis in the lattice Boltzmann equation. In fact, one other basis 
suggests itself naturally. Given a basis of covectors or bras, it is always possible to construct a 
dual basis of vectors or kets. Indeed, our decomposition of the covector space into hydrodynamic 
and kinetic subspaces is naturally mirrored by a corresponding decomposition of the vector space. 
Thus we construct | X a ) where a = 1, . . . , n c and \a v ) where /3 = 1, . . . , n — n c , such that 

(A CT I A r ) = 5 ar (A<j I a v ) = , , 

(a v I A CT ) = {a v I a e ) = 5^. 

Thus, we can expend |N) in terms of these hydrodynamic and kinetic basis kets. 

n c n—n c 

|N> = £ A CT \X a ) + J2 a v K> ■ ( 18 ) 

cr=l T)=\ 

We shall call the set of coefficients X a and the conservation representation. The explicit con- 
struction of this representation is best illustrated by example, and we give several of these in the 
following sections. 

The advantage of the conservation representation for describing collisions is clear: When |N) is 
expanded in terms of hydrodynamic and kinetic kets, collisions may change only the coefficients of 
the latter. From Eqs. (|17|) and (|18|) . we see that the coefficients of the hydrodynamic kets, namely 
the Ao-'s, are the conserved quantities themselves. These are precisely what must remain unchanged 
by a collision. So, in the conservation representation, the collision process is simply an alteration 
of the n — n c coefficients of the kinetic kets, namely the ct^'s. Thus, this representation effectively 
reduces the dimensionality of the space needed to describe the collision process ton- n c . Again, we 
shall illustrate the construction of and transformations between these alternative representations of 
the single-particle distribution function for several examples of entropic lattice Boltzmann models. 

4 Note that we are using the term "hydrodynamic" here to describe degrees of freedom that will result in macro- 
scopic equations of motion, whether or not they are those of a fluid. If mass is the only conserved quantity, a diffusion 
equation generically results; if momentum and/or energy are conserved as well, a set of fluid equations generically 
results. We use the term "hydrodynamic" in either case. This terminology is standard in kinetic theory and statistical 
physics, but often seems strange to hydrodynamicists. 
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2.4 Nonnegativity 

While the conservation representation of |N) is most natural for describing collisions, the original 
representation is more natural for describing the constraint of nonnegativity of the distribution 
function components. In the original representation, we had a restriction to Nj > for j = 1, . . . , n. 
In the conservation representation, these n inequalities transform to a corresponding set of linear 
inequalities on the hydrodynamic and kinetic parameters, A CT and oj„. As is well known, such a set 
of linear inequalities define a convex polytope in the parameter space. This construction, in the full 
n dimensional space of hydrodynamic and kinetic parameters, shall be called the master polytope. 
Specification of the hydrodynamic parameters then defines the cross section of the master polytope 
that bounds the kinetic parameters, a v ; we shall call these cross sections the kinetic polytopes, and 
it is clear that they must also be convex. 

We shall construct the master and kinetic polytopes for simple entropic lattice Boltzmann 
models later in this paper. We shall see that they become very difficult to visualize when the 
number of particle velocities becomes large. This difficulty raises the question of whether or not 
there is a general method to describe the shape of polytopes defined by linear equalities in this 
way. It turns out that such a method is well known in linear programming and optimization theory, 
and is called Fourier- Motzkin elimination [20]. It is constructive in nature, and works for any set 
of inequalities in any number of unknowns. If the inequalities cannot be simultaneously satisfied, 
the method will indicate that. Otherwise, it will yield an ordered sequence of inequalities for each 
variable, the bounds of which depend only on the previously bounded variables in the sequence. If 
it were desired to perform a multiple integral over the polytopic region, for example, the Fourier- 
Motzkin method would provide a systematic procedure for setting up the limits of integration. 

2.5 Collisions 

Once we are able to construct and visualize these polytopes, it is straightforward to describe the 
constraints that conservation imposes on the collision process: Because the propagation process is 
nothing more than a permutation of the values of the Nj(x) on the lattice, it is clear that it will 
maintain nonnegativity. That is, if the Nj(x) were positive prior to the propagation, they will be 
so afterward. This means that the post-propagation state of each site will lie within the allowed 
master polytope. Given such a post-propagation state |N), we transform to the conservation 
represenation. The coefficients X a of the hydrodynamic kets are the conserved quantities and must 
remain unchanged by the collision. Geometrically, these determine a cross section of the master 
polytope within which the state must remain. This cross section is the kinetic polytope of the pre- 
collision state. The essential point is that the post-collision state must also lie within this kinetic 
polytope to preserve nonnegativity. To the extent that the postcollision state is determined by the 
precollision state, this means that: 

Collision Property 1 The collision process is a map from the kinetic polytope into itself. 

We note that this requirement is not without some controversy. It may be argued that lat- 
tice Boltzmann algorithms are ficticious kinetic models from which realistic hydrodynamics are 
emergent. Since the details of the kinetics are ficticious anyway, why not also dispense with the 
requirement that the single-particle distribution function be positive? As long as the conserved 
densities of positive-definite quantities, such as mass and kinetic energy, are positive, why should 
we care if the underlying lattice Boltzmann distribution function is likewise? There are two rea- 
sons: The first is that, even if the system is initialized with nonnegative physical densities, the 
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propagation step may give rise to negative physical densities if negative values of the distribution 
are allowed. To see this, imagine a postcollision state in which all of the neighbors of site x have 
a single negative distribution function component in the direction heading toward x. Even if all 
these neighboring sites had positive mass density, site x will have a negative mass density after one 
propagation step. 

The second reason for demanding nonnegativity of the distribution function is, to some extent, a 
matter of taste. We like to think that the kinetic underpinnings of the lattice Boltzmann algorithm 
are more than just a mathematical trick to yield a desired set of hydrodynamic equations. Though 
there may be no physical system with such a collision operator (not to mention dynamics constrained 
to a lattice), we feel that the more properties of real kinetics that can be maintained, the more 
useful the algorithm is likely to be. This is particularly important for complex fluids, for which 
the form of the hydrodynamic equations is often unknown, and for which we must often appeal to 
some kinetic level of description. Nevertheless, we shall revisit this question in Subsection 15.51 

Collision Property ^ is simple in statement and motivation, but in fact it weeds out many puta- 
tive lattice Boltzmann collision operators, including those most commonly used in computational 
fluid dynamics research today - namely, the overrelaxed lattice BGK operators. For sufficiently 
large overrelaxation parameter (collision frequency), such operators are well known to give rise to 
negative values of the 2V,-(x). This is usually symptomatic of the onset of a numerical instabil- 
ity in the lattice Boltzmann algorithm. We shall discuss such instabilities in more detail later in 
this paper. For the present, we emphasize that we are not claiming that Collision Property ^ will 
eliminate such instabilities. The property does mandate a lower bound of zero on the iVj(x), and 
hence it restricts the manner in which such instabilities might grow and saturate. Nevertheless, it 
is generally still possible for collision operators that obey Collision Property ^ to exhibit numerical 
instability if they lack a Lyapunov functional. 

2.6 Reversibility 

Collision Property^is the minimum requirement that we impose on our lattice Boltzmann collision 
operators. It is possible to define more stringent requirements for them, and we shall continue to 
do exactly that to satisfy various desiderata. For example, we might want to demand that our 
lattice Boltzmann algorithm be reversible. A reversible algorithm could be run backwards in time 
from any final condition by alternately applying the inverse propagation operator, 



followed by an inverse collision operator to recover an initial condition at an earlier time. For 
such an inverse collision operator to exist, the map from the allowed polytope to itself must be 
one-to-one. Since we have already demanded that the map be into, this means that it must also 
be onto, hence: 

Collision Property 2 A collision process is reversible if it is a map from the allowed polytope 
onto itself. 

2.7 Imposition of a Lyapunov Functional 

The criteria that we have set out thusfar assure that the conservation laws - and hence the First 
Law of Thermodynamics - will be obeyed. To ensure stability and thermodynamic consistency, 
however, it is necessary to also incorporate the requirements of the Second Law. 




(19) 
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We suppose that our system has a function H of the state variables Nj that is additive over 
sites x, 

F=£>(N(x)), (20) 
xe£ 

and additive over directions j, 

n 

ft(N(x)) = X;%(JV i ( X )), (21) 
i=i 

where the functions are defined on the nonnegative real numbers. One might wonder if Eq. (|21l) 
could be generalized, but in fact it has been shown to be a necessary condition for an .ff -theorem |21| . 
It is clear from this construction that H is preserved under the propagation operation of Eq. (|13|) . 
If we require that our collisions never decrease the contribution to H from each site - that is, that 
h (N(x)) can only be increased by a collision - then H is a Lyapunov functional for our system, 
and the existence and stability of an equilibrium state is guaranteed. 

We note that h is a function of the Nj for all Nj > 0, and can therefore be thought of as a scalar 
function on the master polytope. Specification of the coefficients X a of the hydrodynamic kets then 
determine the cross section of allowed collision outcomes, or the kinetic polytope, parametrized by 
the coefficients of the kinetic kets a„. Thus, for a given incoming state, h can be thought of as a 
scalar function on the kinetic polytope. We denote this by /i(cO, and we demand that the collision 
process increase this function: 

Collision Property 3 To ensure the existence and stability of an equilibrium state, a collision at 
site x must not decrease the restriction of the function h to the kinetic polytope. 

From a more geometric point of view, note that the (hyper) surfaces of constant ^(o^) provide 
a codimension-one foliation of the kinetic polytope. These (hyper)surfaces degenerate to a point 
where h reaches a maximum. An incoming state lies on such a (hyper)surface. A legitimate collision 
is required to map such an incoming state to an outgoing one that lies inside, or at least on, this 
(hyper)surface. Clearly, the point with maximal h must be mapped to itself. 

Indeed, we can subsume Collision Property ^ into Collision Property 01 by constructing the 
function h so that it takes a minimum value on the boundary of the polytope, and increases to a 
single maximum somewhere inside. The master polytope is defined as the region for which all of 
the iVj(x)'s are nonnegative. The boundary of the master polytope is therefore the place where at 
least one of the Nj(x) , s vanishes. It follows that Ylj Nj(x) is constant - in fact, it is zero - on the 
polytope boundary. More generally, if Cj( x ) f° r J ' = 1> • • • > n are such that 

0(0) =0, (22) 

then YYJ Cj (-^j( x )) a ls° g° es to zero on the polytope boundary. To ensure that there is only one 
maximum inside the convex (master or kinetic) polytope, we also require that the Q be strictly 
increasing, 

Cj(x) > 0, (23) 
for nonnegative x. Indeed, from Eqs. (|22|) and (|23|) it follows that 

Cj(x) > 0, (24) 

for nonnegative x. Thus, the simplest choice would be to make h (N(x)) a function of JX™ Cj (^Y?'( x ))- 
To make this consistent with Eq. 1)2 1|) . however, we see that the functions 6j should be chosen to 
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be the logarithms of the functions Q , so 

n 

MN(x)) = £ln[0(iV,(x))]. (25) 

This goes to negative infinity on the boundaries of the master (and hence kinetic) polytopes, and 
it has a unique maximum in the interior. 



Collision Property 4 A valid functional form for h is given by Eq. \25\) . where the functions Q 
obey Eq. Eq. and Eq. J^p. 

This may well be the most controversial of the four collision properties that we have presented. 
Indeed, it is not necessary to assure that H have a single maximum within the polytope. We shall 
discuss an H function that violates Collision Property 0] in Subsection 15.51 

The equilibrium distribution is the point within the kinetic polytope where h has its maximum 
value. If we denote the kinetic parameters by aj, then we can find this point by demanding that 
the gradient of h vanish, 



m L _ » cm 1 d_N 1 _ » cmi (26) 



If we denote the equilibrium distribution function by Nj q , it follows that the covector with com- 



ponents (Nj^ /Q (NjJ lies in the hydrodynamic subspace, 







(7=1 

The n c coefficients Q a must be chosen so that the values of the conserved quantities are as desired, 

(X T | N cq ) = A r (28) 

for r = 1, . . . , n c . Eqs. (|27|) and (|28|1 can be regarded as n + n c equations for the n + n c unknowns, 
A^ q (for j = 1, . . . , n) and Q a (for a = 1, . . . , n c ). In general, these equations are nonlinear and it 
is not always possible to write the iVj q in closed form. 

The equilibrium distribution is something over which the model builder would like to retain as 
much control as possible, since it is often used to tailor the form of the resultant hydrodynamic equa- 
tions. For example, a judicious choice of the equilibrium distribution function is required to obtain 
Galilean invariant hydrodynamic equations for lattice Boltzmann models of fluids. Customarily, 
lattice Boltzmann model builders have simply dictated the form of the equilibrium distribution 
function. To the extent that this is necessary, however, it is more in keeping with the philosophy 
espoused in this paper to (try to) dictate the form of the Lyapunov functional, rather than that of 
the equilibrium distribution. The challenge is to find a set of functions subject to the require- 
ments of Eqs. (|22j), (|23|) and ((23}, such that Eq. (|27|) yields the desired equilibrium distribution. 
We shall return to this problem in Subsection 15.51 

While this formulation of an entropic principle for lattice Boltzmann models seems reasonable, 
note that it is rather different from the usual one encountered in kinetic theory. The usual choice 
there would be that of Boltzmann, h = j=i Nj In Nj . A moment's examination of Eq. ((2£ 
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indicates that this would correspond to the choice Cj( x ) = % x but this function decreases for 
x £ (0, 1/e) and hence violates Eq. ()23[) . Thus, while it is possible to construct lattice Boltzmann 
models with Maxwell-Boltzmann equilibria the Lyapunov functionals of the models described 
in this paper need to be rather different from those commonly used in kinetic theory. With this 
caveat in mind, we shall henceforth abuse terminology by taking ll H function" and "Lyapunov 
functional" to be synonymous. 

2.8 Collision Operator 

Finally, we turn our attention to the construction of a collision operator that is in keeping with the 
collision properties described above. Obviously, there are many ways to describe maps from the 
kinetic polytopes into (and perhaps onto) themselves, that do not decrease an H function. For the 
purposes of this paper, however, we restrict our attention to the BGK form of collision operator in 
which the outgoing state is a linear combination of the incoming state and the equilbrium, 



where r is called the relaxation time. Usually r is taken to be constant, but more generally it may 
depend on the conserved quantities, and most generally on all of the Nj. 
If we combine Eqs. (|29|) and (j3J), we get the BGK equation, 



From a geometric point of view, this equation tells us to draw a line in the kinetic polytope from 
the position of the incoming state |N) through the equilibrium state |N eq ). The final state is a 
weighted combination of these two states and hence lies on this line. The incoming state is weighted 
by 1 — 1/r, and the equilibrium state is weighted by 1/r. Thus, for r > 1 the outgoing state lies 
somewhere on the segment between the incoming state and the equilibrium. Since both of these 
states are in the kinetic polytope, and since this polytope is convex, the outgoing state must lie 
within it as well. Thus Collision Property^* 13 satisfied. Moreover, since the restriction of h to this 
segment increases as one approaches the equilibrium, Collision Property 01 is also satisfied. 

The instabilities associated with lattice BGK operators arise because practitioners try to over- 
relax them. A large class of lattice BGK models for the Navier-Stokes equations have shear viscosity 
v oc (r — 1/2). In an effort to achieve lower viscosity (and hence higher Reynolds number), prac- 
titioners set r to values between 1/2 and unity. In this situation, the outgoing state "overshoots" 
the equilibrium and lies on the other side of the polytope, opposite the equilibrium state. For 
sufficiently small r, the outgoing state will lie on a contour of h that is lower than that of the 
incoming state, thereby violating Collision Property 01 Still smaller values of r might cause the 
outgoing state to lie outside the kinetic polytope altogether, thereby violating Collision Property^ 
In either case, numerical instability is likely to result. 

2.9 Entropically Stabilized BGK Operators 

A potential solution to this problem is suggested by our geometrical viewpoint. For r = 1 the 
outgoing state is the equilibrium, for which h is maximal. As r is decreased from unity, the final 
value of h decreases from its maximal value. At some value of r less than unity, the outgoing value 




(29) 



T 




(30) 
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of h will be equal to the incoming one. In order to respect Collision Property |3J we must not make 
r any lower than this value, given by the solution to the equation 

h (N') = h (N) , (31) 

where |N') is given by Eq. (|3U|) . This may be regarded as a nonlinear algebraic equation for the 
single scalar unknown r. Indeed, this limitation on r was suggested independently by Karlin et 
al. |22| . and by Chen and Teixeira . It is straightforward to find the solution to this equation 
numerically, since it is easily bounded: We know that the desired solution has an upper bound of 
unity. The lower bound will be that for which the solution leaves the kinetic polytope; this happens 
when one of the iVj's vanish, or equivalently when r is equal to the largest value of 1 — N^/Nj, 
for j G {1, . . . , n}, that lies between zero and one. Given these two bounds on r, the regula falsi 
algorithm will reliably locate the desired solution. Call this solution r*(N). This will generally be 
a function of the incoming state. A useful way to parametrize r is then to write 

r ( N ) = (32) 

where < k < 1 is a constant parameter. The case k — > corresponds to r — > oo so that the 
collision operator vanishes; the case k — > 1 corresponds to r — > r*, which is the largest value 
possible that respects Collision Property [U Thus, the entropically stabilized version of the lattice 
BGK equation is 

|N / ) = |N) + ^-(|N^)-|N)), (33) 

where t* is the solution to Eq. (|31[). and < K < 1. 

Of course, making r a nontrivial function of the incoming state will impact the hydrodynamic 
equations derived from the model. The challenge to the designer of entropic lattice Boltzmann 
models is then to choose the Q very judiciously, so that Eq. P3*|) yields the desired hydrodynamic 
equations, while stability is guaranteed by keeping < k < 1. 

In constructing a lattice Boltzmann model in this fashion, as opposed to using the simpler 
prescription of specifying r, it may be argued that we are relinquishing some control over the 
transport coefficients. After all, if v oc (r — 1/2), it appears that we can specify arbitrarily small 
viscosity by reducing r. In fact, this is not the case since uncontrolled instabilities are known to set 
in for r well above 1/2. The objective of entropic lattice Boltzmann models is to allow the user some 
ability to overrelax the collision process, without sacrificing stability. The price that one may have 
to pay for this stable overrelaxation is living with a bounded transport coefficient. Thus, entropic 
lattice Boltzmann models for fluid flow may be restricted to minimum values of viscosity. As we 
shall show, however, for certain entropically stabilized lattice Boltzmann models, this minimum 
value can actually be zero. This allows for fully explicit, perfectly conservative, absolutely stable 
algorithms at arbitrarily small transport coefficient. 

Finally, we note that this particular prescription is only one way of creating a stable lattice 
Boltzmann algorithm. In fact, any map obeying Collision Properties ^ an d El will work. More 
general mappings of polytopes to themselves can and should be considered, and we leave this to 
future study. 
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3 One-Dimensional Diffusion Model 



Having discussed entropic lattice Boltzmann models in general terms, we now apply the formalism to 
mass diffusion in one dimension. In elementary books on numerical analysis, it is demonstrated that 
the fully explicit finite-difference approximation to the one-dimensional diffusion equation is stable 
only if the Courant condition |24| . At < (Ax) 2 /(2D), is satisfied; note that this places an upper 
bound on the transport coefficient. It is also shown that this condition may be removed by adopting 
an implicit differencing scheme, such as that of Crank and Nicolson [25], or an alternating-direction 
implicit scheme, such as that of Peaceman and Rachford The DuFort-Frankel algorithm |27j is 
fully explicit and unconditionally stable, but it achieves this by a differencing scheme that involves 
three time steps, even though the diffusion equation is first-order in time. 

For the problem of achieving high Reynolds number, one would like the transport coefficient to 
be as small as possible. We shall show that the entropic lattice Boltzmann algorithm provides a 
fully explicit, perfectly conservative, two-time-step algorithm that is absolutely stable for arbitrarily 
small transport coefficient. While this result seems significant in and of itself, part of our purpose is 
pedagogical. In the course of our development of this model, we shall discuss optimal conservation 
representations and present the Fourier-Motzkin method for visualizing the master and kinetic 
polytopes. Nongeneric features of the example are noted, in preparation for more sophisticated 
examples in subsequent sections. 



3.1 Description of Model 

We suppose that we have a regular one-dimensional lattice C whose sites x E C are occupied by 
particles whose velocities may take on one of only n = 3 discrete values, namely —1,0 and +1. We 
abuse notation slightly by letting j take its values from the set of symbols {—,0,+}, so that we 
may write the single-particle distribution as Nj(x). Thus, the state of a given site is captured by 
the ket, 

|N> = N \, (34) 
WW 

where we have suppressed the dependence on the coordinate x for simplicity. 

Diffusion conserves mass, so we suppose that the mass per unit lattice site (mass density in 
lattice units), 

p = iV_ + iV + N + = (p | N) , (35) 
is conserved by the collisions. Here we have introduced the hydrodynamic "bra," 

(p| = ( +1 +1 +1 ) . (36) 

(37) 

This is the one and only hydrodynamic degree of freedom in this example. Because there are a 
total of three degrees of freedom, the other two must be kinetic in nature. To span these kinetic 
degrees of freedom and thereby make the bra basis complete, we introduce the linearly independent 
bras, 

(a\ = ( -1 +2 -1 ) (38) 
</3| = ( +1 -1 ). (39) 
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Next, we form a matrix of the three bras and invert it to get the dual basis of kets, 



( Ip) l«) 1/3) 
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+1 


+i 


+1 


(a) 


- 




-i 


+2 


-1 
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+i 





-1 



-1 




(40) 



From Eq. (|40|) we identify the linearly independent basis kets 




a 



-1 
+2 
-1 



1/3) 




(41) 



The first of these is a hydrodynamic basis ket, while the last two are kinetic basis kets. 

One nongeneric feature should be noted: In this example, we were able to choose (ct\ and 
(J3\ so that each ket is proportional to the transpose of a corresponding bra. Such a choice is 
convenient but unnecessary. The only real requirement in choosing the kinetic bras is that they be 
linearly independent of each other and of the hydrodynamic bras. The next example will illustrate 
a situation in which there is no obvious correspondence between the individual bras and kets. 

We can expand the state ket in this basis as follows 



|N) = p\p) + a\a) + P\P) 



( l-s + M 

1 + a 

\ 1 _ SL — M 
\ L 2 2 



(42) 



where we have defined a = a/p and f3 = (3/ p. The coefficients p, a and 8 (or equivalently p, a and 
(3) constitute the conservation representation. The inverse of this transformation is seen to be 



p = (p | N) = +iV_ + N + N + 
a = (a | N) = -AL + 2N - N + 
P = </3|N) = +N--N+. 



(43) 



3.2 Nonnegativity 

The nonnegativity of the components of the single-particle distribution, |N), places inequality 
constraints on the parameters p, a and /3. For example, it is clear that p > 0. Referring to 
Eq. (|42|) . we see that we must also demand 



o < + M 

2 2 
< 1 + a 

a 38 

< 1 -. 

2 2 



(44) 



It is easy to see that all three of these inequalities may be subsumed by the single statement 

- 1 < a < 2-3|/?| . (45) 
This restricts the kinetic parameters to a triangular region in the (a, 0) plane, as is shown in Fig.^ 
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Polytope for D=1 Diffusion 
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Figure 1: The Kinetic Polytope: Nonnegativity of the distribution function requires that the 
kinetic parameters a and (3 lie in the shaded triangular region. 

Note that the effect of varying p at constant a and (3 is to simply scale the components of the 
distribution function. The triangular region bounding the parameters a and (3 is then independent 
of p. That is, the bounds on the barred kinetic parameters do not depend on the hydrodynamic 
parameter, so there is really no need for distinction between the master and kinetic polytopes. It 
should be noted that this is another nongeneric feature of the present model. While it is always 
possible to scale out by a single nonnegative-definite hydrodynamic parameter, as we have done 
with p in this case, more sophisticated models will have several hydrodynamic parameters. In 
such cases, the shape of the region bounding the kinetic parameters will depend on the remaining 
hydrodynamic parameters. We shall see this more clearly in the example of Section 

3.3 Optimality of Representation 

In this example, there is clearly a great deal of latitude in the choice of the kinetic bras. The only 
requirement that we have placed on these is that they be linearly independent of the hydrodynamic 
bra and of each other. This raises the question of whether there is some optimal choice that might 
be made for these. Of course, this depends entirely on what is meant by "optimal" in this context. 

It has been noted that the original representation of |N) is more natural for expressing the 
constraint of nonnegativity, while the conservation representation is more natural for expressing 
the collision process. For this reason, any computer implementation of entropic lattice Boltzmann 
methods will require frequent transformations between the two representations. This transforma- 
tion is precisely what we have worked out in Eqs. (|42j) and (|43|) above. Thus, one natural figure of 
merit is the number of arithmetic operations required to perform such transformations. 

To investigate this question, rather than requiring that the kinetic bras be given by Eqs. (|38l) 
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and Q39JI . we leave them in the general form 
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(46) 
(47) 



1/3) = 

\ Pi -02 J 
where we have defined the determinant 

A = ori (/? 2 - /%) + «2 (A - /9i) + «3 (A - 02) ■ 




(48) 



(49) 



One way to simplify the transformation process would be to find a representation in which as many 
components as possible of the kinetic bras and kets vanish, without making the transformation 
singular. This means that we want to choose the Oj's and 0^s such that A / 0, while making 
vanish as many as possible of the following twelve quantities: 

ai, a 2 , «3, 
Pi, P2, Ps, 
a 3 — a 2 , a\ — 03, a 2 — a\, 

02- 03, 03 -Pi, P1-P2. 

In this example, it is straightforward to see that there are a number of ways to make six of these 
quantities equal to zero. For example, we could choose 



ai / 0, 



a 2 = 0, 
02^0, 



«3 = 0, 
03 = 0, 



(50) 



resulting in A = a±0 2 7^ and 



a 3 — a 2 = 0, a\ — a 3 7^ 0, a 2 — a\ 7^ 0, 
02-Ps^O, ft -/Ji = 0, Pi -02 7^0. 

This choice also has the added virtue of making the hydrodynamic ket equal to 




(51) 



(52) 



which also has two vanishing components. 

The computation of the |N) components from the parameters p, a and is thus reduced to 



pa 
pP 

p(l-a-P) 
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involving a total of two multiplications and two additions/subtractions; this may be contrasted 
with Eqs. (|42|) which involve six multiplications and five additions/subtractions. Likewise, using 
the choice of Eqs. (|5U|) . the computation of the parameters p, a and (5 from the N is reduced to 

p = N^ + N + N + 
a = AL/p 
= No/p, 



involving two additions/subtractions and two divisions; this may be contrasted to Eqs. (|43|) which 
involve five additions/subtractions, one multiplication, and two divisions. 

Clearly, more sophisticated figures of merit could be devised to optimize the transformations 
used in lattice Boltzmann computations. As we have noted, vanishing components of the kinetic 
bras and kets eliminate the addition/subtraction of terms. Likewise, components equal to ±1 do 
require an addition/subtraction, but not a multiplication, and this fact could be taken into account 
in a more refined figure of merit. The computation of the collision outcome is the principal "inner 
loop" of a lattice Boltzmann computation, insofar as it must be performed at each site of a spatial 
grid at each time step. Such considerations may be especially important for lattice Boltzmann 
models with large numbers of velocities. 

3.4 Fourier-Motzkin Elimination 

For this simple example, we had no difficulty visualizing the master and kinetic polytopes. In 
preparation for the succeeding sections where the task will be substantially more difficult, however, 
we take this opportunity to introduce the Fourier-Motzkin elimination method for this purpose. 
The algorithm consists of the following sequence of steps: 

1. We rewrite the set of inequalities, Eq. (|44[) in matrix format, 

( -\ +1 + 1 \( « A 

+1 0+1 > (53) 

v -s -I + 1 A 1 / 

We have adopted the convention of including constant terms in an extra column of the matrix, 
using the device of appending 1 to the column vector of unknowns. In general, there will be 
m inequalities for n unknowns, and the matrix will be of size m x (n + 1). 

2. We scale each inequality by a positive factor so that the pivot H is either —1, or +1. (Recall 
that scaling by a positive factor preserves the sense of an inequality.) We then reorder the 
inequalities, sorting by their (scaled) pivots so that the zero pivots are last. Beginning with 
Eq. (|HHJ), this yields 

/ -1 +3 +2 \ / a \ 

-1 -3 +2 \\/3 > 0. (54) 
V +1 +1 )\ 1 ) 

Note that there were no zero pivots in this first step. We did, however, reorder the inequalities 
so that those with pivot —1 are together, and preceed that with pivot +1. 

As in discussions of Gaussian elimination, the term "pivot" refers to the first nonzero entry in a row. 
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We now add all pairs of inequalities, such that the first member of the pair has pivot —1 and 
the second member of the pair has pivot +1, and we append the new zero-pivot inequalities 
thus obtained to the system. If there are m_ inequalities with pivot —1, and m + inequalities 
with pivot +1, this results in the addition of m_m+ new inequalities to the system. Since 
there were m — m_ — m + zero-pivot inequalities to begin with, the new system will have a 
total of mo = m — m_ — m + + m_m_|_ zero-pivot inequalities. For our above example, m = 3, 
?n_ = 2 and m + = 1, so we get mj = 2 zero-pivot inequalities in the system, which can now 
be written 

/ -1 +3 +2 \ 
-1 -3 +2 
+1 +1 
+3 +3 
\ -3 +3 / 




> 0. 



(55) 



4. Finally, we recurse by returning to step 2 for the mo x n submatrix obtained by taking only 
the zero-pivot inequalities and neglecting the first unknown (since the zero-pivot inequalities 
do not involve it anyway) . We continue in this fashion until all of the first n elements of the 
rows of the matrix are zero. 

For our example, the recursion step asks us to return to step 2 for the submatrix indicated 
below 
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As before, we begin by normalizing the rows and sorting them by pivot, 
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(57) 



V 

In this example there is now one row of the submatrix with pivot 
append the sum of these to the system to get 

/ 



-1 and one with pivot +1, so we 
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(58) 



V 

At this point, the first n elements of the last row of the matrix are zero, so we stop the recursion 
and examine the last row. It states a true inequality, namely +2 > 0, so we can conclude that the 
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original inequalities are not mutually exclusive; that is, a nonnull polytope of solutions will exist. 
If the n+1 element of the final row had been negative, we would have concluded that no such 
solution was possible. 

Once we have stopped the recursion and established the consistency of the inequalities, we 
can get the final set of inequalities by "back substituting" up the matrix. The last two nontrivial 
inequalities tell us that — 1 < ft < +1, or 

\P\ < 1- (59) 

The first three equations then yield — 1 < a < min (2 - 3/3, 2 + 3/3) which is equivalent to Eq. 
describing the triangular region of Fig. ^ 

More generally, each variable will have a lower bound that is the maximum of all the lower 
bounds determined by the inequalities with pivot +1, and an upper bound that is the minimum of 
all the upper bounds determined by the inequalities with pivot —1. The arguments of the maximum 
and minimum functions will depend only on those variables that have already been bounded. If 
we were using the technique to find limits of integration for a multiple integral, the inequalities at 
the bottom of the Fourier-Motzkin-eliminated matrix would correspond to the outermost integrals. 
We shall revisit this technique in the next section. 



3.5 Collision Operator 

To illustrate the construction of an entropically stable collision operator for this model, we adopt 
the simplest possible H function by taking Cj( x ) = x - Since there is only one hydrodynamic degree 
of freedom, Eqs. (|27j) reduce to 

Nl q = iV cq = Nl q = 1. (60) 
Eq. then tells us that Q = 3/p, so 




(61) 



Comparing this with Eq. (|42|). we see that, within the kinetic polytope, the equilibrium point is 
the origin, a cq = /3 eq = 0. A contour plot of the H function on the kinetic polytope is presented in 

Fig.m 

The lattice BGK equation in the conservation representation, Eq. 1)30(1 . is then equivalent to 
the prescription 

P=P (62) 

(63) 
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where z = 1 — 1/r. These can also be written as the single equation 
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Figure 2: H Function on the Kinetic Polytope: Note that H reaches a maximum when 
a = /3 = 0, and goes to — oo on the boundary of the polytope. 



As promised, the collision is dramatically simplified in this representation. The hydrodynamic pa- 
rameter is unchanged; only the two kinetic parameters are altered by the collision. Transformation 
back to the original representation yields 



|N') 



P 



or equivalently 
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3.6 Entropic Stabilization 

The condition for marginal entropic stabilization, Eq. H31|) . for this model is then 

\nN' + + IniVo + \nN'_ = lniV + + IniVo + lnAL. (67) 

This needs to be solved for the limiting value z = z*, where = 1 — 1/t*. Using Eqs. (|42|) and 
(jHSJ), this becomes 

(i-f +¥)( i+ ^)( i -¥-¥)=( i -l + f)( i+a '( i -f-f)- < 68 > 

This appears to be a cubic equation for z*, but in fact z* = 1 (corresponding to t* — > oo) is clearly 
a root. It is an uninteresting root since it corresponds to no collision. Once it is removed, we are 
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Figure 3: Contours of z*(a,/3) on the Kinetic Polytope: This function approaches —1 when 
a = (3 = 0, is —2 on the midpoints of the sides of the triangle, and is —1/2 at the vertices. 



left with a quadratic for z#, so the relevant solution is easily written in closed form, 



3a 2 - a 3 + 9(3 2 + 9a(3 2 - ^3(a 2 + a 3 + 3/? 2 - 9a/3 2 )(3a 2 - a 3 + 9/3 2 + 9a/3 2 ) 

2(a 3 - 9a^ 2 ) 



(69) 



A contour plot of the function z*(a,j3i) on the kinetic polytope is displayed in Fig. |21 There it 
can be seen that z* has minimum value of —2 at the midpoints of the sides of the triangle, has 
maximum value of —1/2 at the vertices, and approaches —1 at the equilibrium a = f3 = 0. 

The entropically stabilized collision operator is then given by Eq. (|65j). with r — > t*/k, or 
equivalently 

z 1 - k(1 - z*). (70) 
It is useful to discuss three possible choices for k: 

• k — ► means that z — > 1 and r — > oo. This is the uninteresting limit in which the collision 
operator vanishes and hydrodynamic behavior is lost. 

• k = 1/2 means that z = (1 + z*)/2 and r = 2/(1 — z*). This is interesting because near 
equilibrium, where z* ~ —1, it coincides with the prescription r = 1. This is the limit, 
mentioned in Subsection 12.91 wherein the outgoing state is the local equilibrium. This is the 
smallest value of r for which the standard lattice BGK algorithm is guaranteed to be stable. 
(In fairness, the standard algorithm is likely to be stable for smaller r; it's just that this is 
not guaranteed.) 

• k — ► 1 means that z — ► and r — ► r*. This is the largest value of r for which the entropic 
lattice BGK algorithm is guaranteed to be stable. 
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3.7 Attainable Transport Coefficient 



The complicated form of Eq. (|69|) may lead one to believe that the hydrodynamic equation (in 
this case, a diffusion equation) would be very difficult to derive for the entropically stabilized 
collision operator. In fact, this is not the case at all. The Chapman-Enskog analysis that yields 
the hydrodynamic equations requires only that one linearize the collision operator about the local 
equilibrium. Since the local equilibrium has a = (3 = 0, it is clear from Eqs. 1)63|) and (|64j) . that 
only the value of z at equilibrium will enter. Since z* — ► — 1 at the equilibrium, Eq. (|7U|) indicates 
that the transport coefficient will be precisely that obtained by the standard lattice BGK algorithm 
with r = 1/(2k). 

The Chapman-Enskog analysis of the standard lattice BGK algorithm is an elementary exer- 
cise [2Bj- The result for the diffusivity in natural lattice units is 



Hence, by the above argument, the result for the diffusivity of the entropic lattice BGK algorithm 
is 



From this, we can clearly see the benefit of the entropic algorithm. The standard lattice BGK 
algorithm is not guaranteed stable for r < 1 or k > 1/2, as discussed above. From Eq. ()71j) or 
(17211 . respectively, we see that this corresponds to a minimum diffusivity of 1/6. By contrast, the 
entropic lattice BGK algorithm loses stability only when k = 1, corresponding to D = 0. 

It is remarkable that we have found a stable, conservative, explicit algorithm that seems to 
allow the transport coefficient to become arbitrarily small. After all, the analysis leading to the 
shear viscosity of a fluid model is not very different from that leading to the diffusivity above, and 
arbitrarily small shear viscosity would allow for arbitrarily large Reynolds number. In the field 
of computational fluid dynamics, however, results that seem too good to be true usually are just 
that. For one thing, stability does not imply accuracy. A perfectly stable algorithm is not terribly 
useful if it does not converge to the correct answer. Flows at higher Reynolds number involve ever 
smaller eddies, and at some point the lattice spacing becomes insufficient to resolve these. While 
it is comforting to have an algorithm that does not lose stability in this situation, it is probably a 
mistake to attach much physical significance to its results. 

Another problem is that the time required for the system to come to equilibrium goes to infinity 
as k — ► 1. In fact, when k = 1 the entropic lattice BGK collision operator is its own inverse: If 
incoming state |N) yields outgoing state |N'), then incoming state |N') would yield outgoing state 
|N). This means that if the entire lattice were initialized slightly away from equilibrium with no 
spatial gradients whatsoever, it would simply thrash back and forth between two states, without 
ever converging to the desired equilibrium. When the time required to achieve local equilibrium 
exceeds time scales of interest, hydrodynamic behavior breaks down. 

To illustrate this problem, we simulated the entropically stabilized diffusion model with the 
initial density profile 



on a periodic lattice. We used lattice size L = 32, and wavenumber t = 3 in all the simulations 
reported below. We initialized the state of each site in a local equilibrium iV_ = Nq = N + = p/3, 
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Figure 4: Decay of sinusoidal density profile for k = 0.9 (left) and k = 0.99 (right). 
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Figure 5: Decay of sinusoidal density profile for k = 0.999 



without the Chapman-Enskog correction due to gradient. Now Eq. (|73j) is a solution of the diffusion 
equation, if p\ decays in time as exp(— jt), where 



D 



(74) 



So we fit our results to the functional form of Eq. (J73J), measured pi(t), and did a least-squares 
fit of its logarithm to a linear function of time in order to determine 7. We then used Eq. (|74|) to 
get the diffusion coefficient, and we compared this to the theoretical value provided by Eq. ((721) for 
several different values of k, approaching unity. 

Fig. 0] shows the measured value of pi(t) for k = 0.9 and k = 0.99. As can be seen, there is 
an initial transient, due to the inadequacy of the form used for the local equilibrium iV_ = Nq = 
N + = p/3 in the presence of a spatial gradient. The left-hand side of Fig. shows the measured 
value of pi(t) for k = 0.999, and the right-hand side is an enlargement of the transient region. 
Fig. IHlshows the same things for k = 0.9999. Note that the transient period lengthens considerably 
as k nears unity. We could have substantially reduced these transient periods had we used the 
Chapman-Enskog correction to the local equilibrium; instead, we simply waited for the transient 
to die away before measuring the decay constant 7. 

The results for the diffusivity are displayed in Tabled Note that we find reasonable agreement 
until we get to the cases k = 0.9999 and k = 0.99999. To see what is going wrong for these cases, 
the upper-left-hand corner of Fig. shows the measured value of pi(t) for k = 0.99999, and the 
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Figure 6: Decay of sinusoidal density profile for k = 0.9999 
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2.28751 x 10- 


-6 



Table 1: Theoretical and measured values of the diffusivity for various values of ac. 



upper-right-hand corner shows the same plot with a reduced range for the ordinate. In addition 
to the oddly shaped envelope of the initial transient period, similar to that seen in Fig. |HJ and 
enlarged in the lower part of this figure, we note that the plot of pi(t) in the upper right seems 
to be increasing in thickness, even after the initial transient. Upon further enlargement, shown on 
the left-hand side of Fig. |H1 we see that this thickness is in fact a subtle transient of extremely long 
duration. The right-hand side of the figure shows this same transient at the longest time for which 
the simulation in Fig. [7| was run. Though reduced in magnitude (the scales of the ordinates of both 
graphs in Fig. [HI are equal), the transient is still present, even though most of the signal itself has 
diffused away. We believe that this transient is causing the anomaly in the measured diffusivity. A 
much longer simulation would be required to test this assertion. 

To better understand this loss of separation between the kinetic and hydrodynamic time scales, 
note that the time required for the signal to decay away is 

11k 

Signal OC - OC — OC . (75) 

7 D 1 — ac 

The time required for the transient to decay to 1/e of its initial value, on the other hand, obeys 

K r transient = _ (76) 

e 

whence 

^transient K , • (77) 

In At 

For k = 1 — e, where e is small, the characteristic times in Eqs. ((75)1 and (|77|) both go like 1/e. 
Thus, the transients begin to linger for hydrodynamic time scales, and it becomes impossible to 
separate the kinetic behavior from the hydrodynamic behavior. 
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Figure 7: Decay of sinusoidal density profile for k = 0.99999 
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Figure 8: Decay of sinusoidal density profile for k = 0.99999 



29 



When the time required for the decay of transients is comparable to hydrodynamic time scales 
of interest, the usefulness of the simulation is questionable. Thus, while Table ^ indicates good 
agreement with theory, and the transport coefficient really does tend to zero without loss of stability, 
one must be exceedingly careful in the preparation of the initial condition to exploit this feature of 
entropic lattice Boltzmann models. In particular, transient behavior could be dramatically reduced 
by including the Chapman-Enskog correction to first-order in the gradient. Indeed, if this problem 
turns out to be the only obstacle to stable, conservative, explicit algorithms with arbitrarily small 
transport coefficients, higher-order gradient corrections may also be worth investigating. 



30 



4 One-Dimensional Compressible Fluid Model 



In this section, we apply the entropic lattice Boltzmann method to a simple five- velocity model of 
fluid dynamics in one dimension, first considered by Renda et al. in 1997 |T3|. We shall find that 
the geometric picture is much richer than that for the diffusion model. The master polytope for 
this model is four-dimensional, and we shall show that the Fourier-Motzkin algorithm is very useful 
for describing it. 

4.1 Description of Model 

As a second example, we consider a lattice Boltzmann model for a one-dimensional compressible 
fluid with conserved mass, momentum and energy, first studied by Renda et al. [Tl|- The velocity 
space of this model consists of five discrete values of velocity, namely 0, ±1, and ±2. The single- 
particle distribution at site x with velocity j G {—2, —1,0, +1, +2} is denoted by Nj. As in the last 
example, the state vector |N) at a given site x will be denoted by a ket, 



where we have suppressed the dependence on x for simplicity. 

A compressible fluid conserves mass, momentum and kinetic energy, so we suppose that the 
corresponding densities 



N) 



I N - 2 \ 

JV_! 

N 
N+i 
\ N +2 J 



(78) 



P 



E N i = (p I N > 



(79) 



i=-2 

+2 



TT 



E 3 Nj = (tt I N) 



(80) 



£ 




(81) 



are conserved by the collisions. Here we have introduced the bras 
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or explicitly 
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+1 +1 +4 



(85) 



These are the three hydro dynamic degrees of freedom in this example. Because there are a total 
of five degrees of freedom, the other two must be kinetic in nature. To span these kinetic degrees 
of freedom and thereby make the bra basis complete, we introduce the linearly independent bras, 

(a\ = - ( +1 0+1 ) (86) 
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(87) 



Next, we form a matrix of the five bras and invert it to get the dual basis of kets, 
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From Eq. (|89j) we identify the linearly independent basis kets, consisting of the hydrodynamic basis 
kets 

/ \ ( -1 \ / +1 \ 
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V +i J 



(90) 



and the kinetic basis kets 
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(91) 



Unlike the previous example, there is no obvious correspondence between the individual bra and ket 
basis elements. This is because there is no natural notion of a metric in this space. Nevertheless, 
as noted in Section |2 it is always possible to construct an entire ket basis from an entire bra basis, 
and that is what we have done here. 

We can then expand the state ket in this basis as follows 
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P\P) 
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( \ {£-Tt-a + 20) \ 
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(92) 



where we have continued to use an overbar to denote conserved quantities per unit mass; for 
example, 7f = 7r/ p is the hydrodynamic velocity. As with the last example, this representation of 
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the distribution function |N) has the virtue of making the conserved quantities manifest. In this 
case, we have the inverse transformation 
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4.2 Nonnegativity 

We now consider the shape of the master and kinetic polytopes for the compressible fluid model 
described in the last subsection. Here we shall find a much richer structure than the simple triangu- 
lar region that we found for the diffusive example in the last section above. Referring to Eq. (|92|) . 
we see that that nonnegativity of the components of the single-particle distribution function is 
guaranteed by the set of inequalities p > and 

< e -it -a + 2(3 
< a-P 

< 2-e-3a (94) 

< a + (3 

< e + it -a -2/3. 

These define the master polytope in the four dimensional (tt, e, a, (3) space. If we fix the hydro- 
dynamic parameters, tt and e, the above still constitute a set of linear inequalities specifying the 
corresponding kinetic polytope in the two dimensional (a, (3) space. Geometrically, these are the 
intersections of the master polytope with the (hyper) planes of fixed hydrodynamic parameters 7f 
and e. Thus, though the mass density p still scales out of the distribution function, the shape of the 
kinetic polytopes in the (a, (3) plane will depend on the hydrodynamic velocity 7f and the kinetic 
energy per unit mass e. 

It is not particularly easy to visualize the shape of the kinetic polytope in the (a, (3) plane 
for given hydrodynamic parameters 7f and e. In fact, the task is tedious enough that we have 
relegated it to Appendix El which takes the reader on a detailed tour of the four-dimensional 
master polytope and its kinetic polytope cross sections. To summarize the conclusions of that 
Appendix, the projection of the master polytope on the (jr,s) plane is illustrated in Fig- El Values 
of (tt, e) outside the shaded regions are not possible. The shape of the kinetic polytope is then a 
triangle when (tt , e) is in the most lightly shaded regions of Fig. ^3 a quadrilateral when (tt, e) is 
in the intermediately shaded regions of Fig. |H1 and a pentagon when (tt, e) is in the most darkly 
shaded region of Fig- EI Illustrations and detailed descriptions of the kinetic polytopes for particular 
values of (it, e) in each of the seven distinct regions with it > are provided in the seven figures 
of Appendix El the chart below Fig. El indicates which figure corresponds to each of its regions. 
Finally, since the inequalities are invariant under (it, (3) — > (—it, —(3), the kinetic polytopes for it < 
are obtained from those for it > by simply reflecting in (3. 
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Figure 9: Bounds on tt and e: The regions with the lightest shading correspond to triangular 
regions in the (a, (3) plane, those with intermediate shading correspond to quadrilateral regions in 
the (a, (3) plane, and that with the darkest shading corresponds to pentagonal regions in the (a, (3) 
plane. 

4.3 Fourier-Motzkin Analysis for the ID Compressible Fluid Model 

To compute the master polytope for the compressible fluid model, we cast the inequalities of 
Eqs. (|94|) into matrix format as follows, 
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Here we have put the hydrodynamic parameters tt and e after the kinetic parameters a and (3 in 
the column vector of unknowns, for reasons that will become apparent in a moment. Performing 
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the Fourier-Motzkin algorithm on this set of inequalities, we finally arrive at the 23 inequalities, 
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(96) 



The last of these is a consistency condition which is seen to be satisfied. The five inequalities prior 
to that yield 

l 



max 



, -1,-4| <e< 2, 



or 



< e < 2, (97) 
which explains the extent of Fig. |U] along the ordinate. The six inequalities prior to that yield 
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(98) 



which explains the extent of Fig. along the abscissa. Thus, we have bounded the projection of 
the master polytope in the hydrodynamic parameter space. 

Given hydrodynamic parameters in the allowed region thus determined, the first 11 inequalities 
of Eq. (|96|) then specify the kinetic polytopes. (This is why we put the hydrodynamic parameters 
after the kinetic ones in the column vector of unknowns.) Continuing our back substitution, the 
bounds on (3 are 
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(99) 
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Finally, the bounds on a are 

max [p, -ji) < a < min (+2/3 -n + s, -2(3 + vf + e) . (100) 

Eqs. and (|100j) specify the shape of the kinetic polytope in the (a, (3) plane, given the hydrody- 
namic parameters n and e. In fact, they provide a more succinct way of determining and expressing 
the bounds of both the master polytope and the kinetic polytopes than the analysis presented in 
Appendix 1X1 
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5 Two and Three-Dimensional Fluid Dynamics 



In this section, we apply our methodology to much richer models, namely those most commonly 
used for two and three-dimensional Navier-Stokes hydrodynamics. Thought the Fourier-Motzkin 
method cannot easily be applied to the latter, we are able to present a valid collision representation 
for both examples. We indicate how a computer simulation for such models might be implemented. 



5.1 FHP Model 

The diffusion and compressible fluid models that we considered above had three and five particle 
velocities, respectively. Such models are of academic and pedagogical interest, but not particularly 
useful as simulational tools. The smallest model that has been used for serious computational fluid 
dynamics calculations in two dimensions is called the FHP model. This is a six-velocity model on 
a triangular grid which has been shown to yield isotropic, incompressible Navier-Stokes flow in the 
appropriate scaling limit P3 12] • The six velocities are given by 
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for < j < 5. The model conserves mass with density 
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and momentum with density 



so the hydrodynamic bras are 
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The kets can then be found as with the preceeding examples; the hydrodynamic kets are 
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and the kinetic kets are 
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The general representation of the state of the FHP model is then 

|N> = p\p) + n x \iv x ) +K y \Tr y ) + a\a) + fi\(3) + 

/ l + 2n x + a-2P \ 

1 + 7f x + V^Tty - a + (3 + V3j 
1 — 7f x + \/3% + a + /?_- \Z3~7 
1-2^-0-2/3 

1 — 7f x — a/3% + a + P + a/37 
y 1 + tt x - \/3Tt y - a + (3 - \/37 y 

In addition, the factors of y/3 could be absorbed into Tt y and 7 to further simplify the final expression. 
5.2 FCHC Model 

Early attempts to generalize the FHP model to obtain a three-dimensional lattice model of the 
Navier-Stokes equations were thwarted by the observation that no regular three-dimensional lattice 
would yield the required isotropy - as the triangular lattice did in two dimensions. The problem 
was solved by noting that a lattice with the required isotropy existed in four dimensions, and its 
three-dimensional projection could be used to construct an isotropic model 0. The four dimen- 
sional lattice is called the Face-Centered Hypercubic (FCHC) lattice. It is a self-dual lattice with 
coordination number 24. It can be defined as all sites (i,j,k,l) on a Cartesian lattice such that 
i + j + k + 1 is even. The 24 lattice vectors - or, equivalently, the 24 lattice sites that neighbor the 
origin (0,0,0,0) - are enumerated in Tabled 

To project these lattice vectors to three-dimensional space, we simply ignore the ficticious I 
coordinate. When we do this, we note that 6 of the resulting three-dimensional lattice vectors will 
have two preimages from the original set of 24 FCHC lattice vectors. Rather than maintain these 
as separate three-dimensional vectors (as is done in lattice-gas studies), it is possible to combine 
them in lattice Boltzmann studies, simply weighting contributions from those directions by a factor 
of two. In fact, this weighting factor can be thought of as a direction-dependent particle mass rrij, 
so that the three-dimensional projection of the FCHC model is as presented in Table 01 

The model is required to conserve mass and three components of momentum, so that the four 
hydrodynamic bras are given by 

(114) 
(115) 
(116) 
(117) 

These are shown in the first four rows of Table 0] in Appendix El followed by a choice of twenty 
linearly independent kinetic bras. The corresponding kets are shown in Table El in Appendix El 
and these can be used to construct a general analytic form for the state ket, |N). 
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Table 2: Lattice vectors of the FCHC model. 
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Table 3: Lattice vectors and masses of the three-dimensional projection of the FCHC model. 



5.3 Attainable Transport Coefficient 



As with the one-dimensional diffusion model, we expect that the entropic lattice Boltzmann model 
will allow us to reduce the transport coefficient of the model. For fluids, this would mean smaller 
shear viscosity and hence higher Reynolds number. To see if this is possible, we consider the FHP 
model in the incompressible limit. In the incompressible limit, the Mach number scales like the 
Knudsen number, so the Chapman-Enskog method involves only the equilibrium distribution at 
zero momentum. Taking Cj( x ) = x t h is straightforward to extremize h to find the equilibrium 
distribution. This may be done perturbatively in Mach (equivalently, Knudsen) number, and the 
result to tenth order is 
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(118) 

(119) 
(120) 



Thus the equilibrium values of (3 and 7 are of second order, while that of a is actually third order. 
For the computation of viscosity in the incompressible limit, the Chapman-Enskog analysis requires 
the linearized collision operator only to zeroth order in the Mach number. Thus, for the purposes 
of working out the attainable viscosity, we may simply set the equilibrium values of all three kinetic 
parameters to zero. It follows that the outgoing kinetic parameters are equal to the incoming ones 
multiplied by z = 1 — 1/r. 
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We can now construct an entropic lattice Boltzmann model. Using Eq. ()113l) for guidance, we 
see that the analog of Eq. (|68|) for is 

(l + 2tt x + a - 20) (l +7r x + V3tt v -a + P + V3 7 ) (l - 7f a + \fM v + a + fi- V3j) 
(l - 2it x - a - 20) (l-it x - \/3n y + a + j3 + v 7 ^) (l + 7T* - V3tv v - a + /3 - y/3j) = 

(l + 2n x + z*a - 2z*0) (l + n x + \/3iiy — z*a + z»(3 + \/3z*j) (l — 7T X + V^^y + z*a + z*f3 - y/iz*^) 

(l - 2n x - z,a - 2z,0) (l - n x - s/3ity + z t a + z*/3 + \/3z^) (l + ty x - \/3iTy - z,a + z*0 - V3z t j) , (121) 

where it bears repeating that, just for the purposes of this computation of viscosity, we have 
neglected corrections to the equilibrium of order Mach number squared H. This is a sixth-order 
equation for z*; since z* = 1 is clearly a root, it can be reduced to fifth order. In spite of the fact 
that fifth order equations generally do not admit solutions in radicals, this one does happen to do 
so. The result, however, is complicated and uninspiring, and so we omit it here. The interested 
reader can simply type the above equation into, e.g., Mathematica® , and use the Solve utility to 
see the result. The important point is that the relevant root for approaches — 1 as a, (3, 7 — > 0. 

The rest of the analysis is exactly as it was for the one-dimensional diffusion model. We set 
z = 1 — k(1 — z*), where < k < 1. Since the usual lattice Boltzmann version of the FHP model 
has the shear viscosity proportional to r — 1/2, the entropic version will have it proportional to 
1/k — 1 which goes to zero as k approaches unity. Moreover, absolute stability is maintained as 
this limit is approached, even though the algorithm is fully explicit and perfectly conservative. As 
with the diffusion model, we are prevented from achieving arbitrarily small transport coefficient 
only by considerations of accuracy (resolution of the smallest eddies), and of the time scale for the 
approach to equilibrium. 

5.4 Computer Implementation 

We pause to consider the computer implementation of the entropic lattice Boltzmann algorithm 
for the FHP fluid. The propagation step can be handled in precisely the same fashion that it is for 
conventional lattice Boltzmann simulations. The principal difference lies in the implementation of 
the collision operator, and that is what we describe here. Unlike the analysis in the last subsection, 
a complete computer implementation needs to solve for the equilibrium exactly - not just to within 
the Mach number squared. 

Given the six quantities, Nj for j = 0, . . . , 5, entering a site, we first contract them with the 
bras, Eqs. (|lU4j) through H1U9|) . to obtain the conservation representation, namely p, ir x , ir y , a, (3 
and 7. Given tt x and tt v , it is then necessary to solve for the equilibrium values of a, (3 and 7. 
We can use Eqs. ()118j) through (|12U|) to get excellent approximations for these, and then use a 
Newton-Raphson iteration to get the correct values to machine precision. Even with an excellent 
initial guess, there is some potential for loss of convergence in any Newton-Raphson iteration, and 
so it might be worthwhile to investigate better ways of solving for these values. Assume that we 
can do this, and denote the equilibrium values with an "eq" superscript. 

Given the equilibrium values, we set 



Note that these terms of order Mach number squared are necessary to derive the remainder of the Navier-Stokes 
equations. 
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and find the value of r that solves the nonlinear equation /i(N') = /i(N). This equation will look 
much like Eq. Ijl21j) . except that <5 eq , /? eq and 7 eq will enter, corresponding to corrections of order 
Mach number squared. Since this nonlinear equation has only one scalar unknown, namely r, it 
can be solved using an iteration method that is guaranteed to converge, such as regula falsi. We 
know that r = 1 is one bound on the solution; the other bound can be taken as the point at which 
the line of outgoing states intersects the boundary of the kinetic polytope. This means finding the 
largest value of r between zero and one for which N'- = for some j G {0, . . . , 5}. Using these two 
bounds, we iterate the regula falsi algorithm to obtain the solution r*. We then set r = t*/k, and 
use Eqs. (|122|) through (|124j) to get the outgoing values of the kinetic parameters. We finish the 
collision procedure by using Eq. Q113[) to get the outgoing state in the original representation. 

5.5 Galilean Invar iance 

As noted in the Introduction, practitioners of lattice BGK models have long known that it is 
necessary to tailor the equilibrium distribution function in a certain way in order to maintain 
Galilean invariance. Recall that the zeroth moment (with respect to the velocity vector) of the 
distribution function is determined by the mass density. Likewise, the first moment is determined 
by the momentum density. For a lattice BGK model of an ideal gas |29| I30j . the trace of the second 
moment is determined by the kinetic energy density. To obtain the correct, Galilean-invariant form 
of the compressible Navier-Stokes equations, it is also necessary |31j to mandate the rest of the 
second moment, the third moment, and the trace of the fourth moment. 

While it is possible to construct equilibrium distributions with these moments for various lat- 
tices, no lattice model has ever been found for which these equilibria are guaranteed to be nonneg- 
ative. This means that the equilibrium toward which we would like to relax in order to maintain 
Galilean invariance may itself lie outside of the master polytope! There is thus no way that such 
an equilibrium could be the maximum of an entropy function that goes to — oo on the polytope 
boundary, such as do the ones that we have been considering. 

In such a situation, we may ask whether it is possible to jettison the requirement of nonnegativ- 
ity, while still demanding absolute stability. In fact, our formalism suggests an interesting way of 
doing just that. If we no longer demand that the function h go to negative infinity on the polytope 
boundary, then we need no longer demand Eq. (|22|). though Eq. I|23ft is still useful to ensure that 
the equilibrium is unique, and Eq. (|2*4*|) is useful to keep h real- valued in light of Eq. (|25|) . Thus, 
we are free to make the choice Cj( x ) is ex P[ rz;2 /(25j)], where the g^s are positive constants. Indeed, 
this choice is similar to one made by Karlin et al. [H2]- Eq. (|27|1 is then the set of n equations, 



and Eqs. (|28|) give n c additional requirements. One may then try to solve these n + n c equations 
for the n + n c unknowns, Q a and gj, subject to the requirement that gj > 0. In fact, not all of 
these unknowns are independent; the scaling gj — > gj/x and Q a — ► Q a %, where x is arbitrary, leaves 
Eq. (|125|) invariant, so there are more requirements than equations, and hence not all equilibria are 
derivable from an h function of this form. Nevertheless, many interesting equilibria are derivable 




(125) 
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in this way, and this approach leads immediately to the construction of an entropically stabilized 
algorithm for such equilibria. Investigation of this approach is work in progress [33] • Its successful 
resolution may yield an absolutely stable entropic lattice Boltzmann algorithm for compressible 
flow, albeit one that allows for negative Nj. 
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6 Conclusions 



In summary, we have presented a general methodology for constructing lattice Boltzmann models 
of fluid dynamics that conserve a given set of moments of the distribution function. The method 
guarantees the nonnegativity of the distribution function, and allows for convenient visualization of 
the dynamics by construction of the polytope of allowed states using Fourier-Motzkin elimination. 
We have shown how such models can be endowed with a Lyapunov functional, resulting in uncondi- 
tional numerical stability, and we presented various choices for the H function for this purpose. We 
also described the computer implementation of such models in some detail. For lattice Boltzmann 
models of diffusion and of incompressible fluid flow, we were able to present fully explicit, perfectly 
conservative, absolutely stable algorithms that appear to work for arbitrarily small transport coef- 
ficient. Indeed, we showed that the limitations on attainable transport coefficients for these models 
arise from considerations of accuracy, rather than stability. 

There are numerous avenues for extension of the current work. We close by presenting a list of 
these. We hope that the richness and promise of this style of model will inspire others to take up 
some of these questions. 

(i) We have not discussed the question of boundary conditions in this paper. For a viscous fluid, 
for example, collisions at a solid boundary wall must conserve mass, but not momentum or 
energy. The latter two quantities can enter or leave the domain through the wall, resulting 
in drag, lift, surface heating, etc. This results in very different polytope structure at the 
boundaries than in the interior. In addition, there is nothing preventing entropy from entering 
or leaving through the wall, so the entire question of how to construct h functions at the 
boundary sites needs to be revisited. 

(ii) While we touched on the question of optimality of conservation representations, clearly much 
more work could be done along these lines. In addition to making it computationally easier 
to transform back and forth between the original and conservation representations, one might 
like to try to optimize the form of the entropically stabilized collision operator. 

(iii) We mentioned reversible models, but did not discuss or study them in any detail. Our 
entropically stabilized models are reversible in the limit as k — > 1, but this is precisely the 
limit in which they begin to fail to simulate the parabolic equations of interest, namely the 
diffusion and Navier-Stokes equations. It would be interesting to know if such reversible 
models are capable of simuating time-symmetric partial differential equations, such as Euler's 
equations of inviscid fluid dynamics. 

(iv) The nature of the anomalies encountered as k — > 1 remains to be completely elucidated. For 
example, the lower right plot in Fig. [7| appears to indicate that the approach to equilibrium is 
extremely slow for this model. Does it eventually reach an exponentially decaying equilibrium 
state for which the diffusion coefficient matches that of the theory? Does roundoff error play a 
role in this regime? These open questions will require substantially more numerical simulation 
than that presented here. 

(v) Given that the anomalies as k — > 1 are due to the slowness of the approach to the equilibrium 
state, much effort could be expended in trying to start the simulation in as close to a global 
equilibrium state as possible. This means taking spatial gradient corrections into account. At 
a minimum, the Chapman-Enskog corrections - accurate to first-order in the spatial gradient 
- could be incorporated with little effort. 
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(vi) The question of how to deal with interaction potentials in this formalism could be addressed. 
At the level of mean-field theory, this would require the introduction of conserved quantities 
that are quadratic in the distribution function. The allowed regions would then no longer 
be polytopes, convex, or even simply connected. Moreover, multiple local equilibria within 
the allowed region may be expected. For near-equilibrium dynamics, this might give rise to a 
Ginzburg-Landau model with multiple free-energy minima. This is thus a difficult extension 
of the current theory, but one for which the rewards would be great. 

(vii) One could try to construct an adaptive mesh refinement (AMR) version of this algorithm 
for which the nonnegativity of the physical densities is the refinement criterion. Using such 
an approach, it might be possible to relax the criterion of nonnegativity of the distribution 
function. That is, we could carry negative distribution function values, but insist that all 
physical densities be positive. Should the propagation step threaten to create a negative 
density, we would try to refine the lattice locally until that is no longer the case. 

(viii) Finally, it would be interesting to consider the potential utility of this algorithm for the con- 
struction of "eddy-viscosity" subgrid models of turbulence. As noted at several points in this 
paper, while entropic lattice Boltzmann algorithms may be stable for arbitrarily small trans- 
port coefficient, they may begin to lose accuracy. The smallest eddy sizes in three dimensional 
Navier-Stokes turbulence, for example, scale as Re -3 / 4 , where the Reynolds number Re goes 
inversely with viscosity. Thus, at sufficiently small viscosity, the eddies will be smaller than any 
fixed grid spacing. If one's goal is an ab initio simulation of the Navier-Stokes equations, this 
circumstance is grounds for rejecting the entropic lattice Boltzmann solution as unphysical. 
On the other hand, stable, coarse-grained, turbulent flow must obey the same conservation 
laws as the underlying Navier-Stokes equations. Since the entropic lattice Boltzmann colli- 
sion operators are the most general that respect these conservation laws, it is interesting to 
speculate as to whether some physical significance might yet be attached to their solutions, 
even when the smallest eddies are not resolved by the lattice. 
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A A Tour of the Master Polytope for the Compressible-Fluid 
Model 



The five inequalities of Eq. (|94|) can be recast as 

< \0\ < a < min {^^,e-\2(3 - vf^ . 

This implies 

e > a+ |2/3-vf| > a > \(3\ > 

and 

e < 2 -3a < 2-3|/3| < 2, 

and hence a > and e G [0,2]. Moreover, we note that all of the inequalities 
the transformation 

7f — > — it 

- -A 

so that we may restrict our attention to it > without loss of generality. 

Next note that the second argument of the min function in Eq. (|126[) will be smaller than the 
first if < e < (2 — e)/3, or < e < 1/2. In this case, for sufficiently small it, the boundary of the 
polytope in the {a, pi) plane will be a quadrilateral, such as that in Fig. 1101 whose lower bound in 
a (shown in black) is \P\ and whose upper bound in a (shown in gray) is e — 1 2/5 — vr | . The allowed 
region of the (a, (3) plane is the shaded area in between these bounds. The bottom vertex (a) of 
this quadrilateral is at the origin, (a,P) = (0,0), and the upper vertex (c) is at (a,P) = (e, it/2). 
The right vertex (b) is then at (3 = e— (2(3 — it) or (a,P) = ((e + 7f)/3, (e + 7f)/3), and the left vertex 
(d) is then at — (3 = e + (2/3 — tt ) or (a, p) = ((e — 7f)/3, (— e + 7f)/3). These results are summarized 
in the table included in Fig. ITUl 

We now ask for what range of e and w the above-described quadrilateral boundary in the (a, P) 
plane is valid. We note that vertices (a) and (b) will be degenerate when e = — it, and that vertices 
(a) and (d) will be degenerate when e = +tt. Also, it is easy to see that one or the other of 
these two degeneracies will occur before any degeneracy involving vertex (c). It follows that the 
above-described quadrilateral boundary in the (a,/3) plane is valid only for |7f| < e < 1/2. This 
region of the (it, e) plane is the shaded triangle AIH in Fig. |H1 

Beginning in triangle AIH, boundary AI is encountered when e = it so that vertices (a) and 
(d) are degenerate. If we cross this boundary the allowed region in the (a, P) plane is no longer 
a quadrilateral, but rather becomes a triangle, such as that in Fig. 1111 whose upper bound in a 
(shown in gray) is s — \2/3 — it\ and whose lower bound in a (shown in black) is (3. The bottom 
vertex (e) of this triangle is at a = (3 = e + (2(3 — it) or (a, (3) = (it — e, it — e); this vertex replaces 
the degenerate vertices (a) and (d) of the above-described quadrilateral. The expressions for the 
coordinates of the triangle's upper vertex (c) and its right vertex (b) are identical to those of the 
corresponding vertices of the quadrilateral. These results are summarized in the table included in 
Fig. El 

Again, we ask for what range of e and it the above-described triangular boundary in the (a, (3) 
plane is valid. We note that vertices (b) and (c) will be degenerate when (e + 7f )/3 = e or e = it/2; 

7 Throughout this appendix, we denote vertices in parentheses, so as not to confuse them with other variables. 



(126) 

(127) 
(128) 

are invariant under 
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for smaller values of e the set of allowed points in the (a, (3) plane is empty. We also note that, just 
as for the quadrilateral of Fig. I1U1 the first argument of the min function in Eq. ()126|) will become 
the determining factor when e > 1/2, also invalidating the above argument. Thus, the allowed 
region in the (a, p) plane will be a triangle with vertices described in Fig. only if (e, ff) is in the 
shaded triangle ABI in Fig- EI Since the region AEH is the image of ABI under the map n — ► — 7f, 
it follows that the allowed region of the (a, (3) plane is also triangular for (e, 7f) G AAEH, with 
vertices identical to those in Fig. ^Jwith (3 — > —(3. 

We have now described all of Fig. below the line s = 1/2 (that is, below line EB). We 
next consider the situation for e > 1/2. This means that we have to start taking into account the 
inequality 

a < ^ (129) 

of Eq. (|126|) . This is a constant upper bound on a which will become less than the a coordinate 
of vertex (c) in Figs. El an d ^2 when e > 1/2. This results in a horizontal truncation of the 
quadrilateral of Fig. El so that it becomes a pentagon, and of the triangle of Fig. ^2 so that it 
becomes a quadrilateral. These are shown in Figs. ll2l and ll31 respectively. The top vertices (/) and 
(<?) in both of these figures are at the upper bound a = (2 — e)/3, and (2 — e)/3 = e =F (2/3 — v) or 
(3 = (3v ±4e^2)/6, respectively; these vertices replace vertex (c) in Figs-Eland^D The expressions 
for the coordinates of vertices (a), (6), (c£) and (e) are identical to those derived previously. The 
vertex coordinates are given in the tables in their corresponding figures. 

Yet again, we ask for what range of e and tt the boundaries pictured in Figs. El and El are 
valid. As e increases, the upper bound on a decreases until vertices (b) and (/) coincide. This 
happens when (e + vf)/3 = (2 — e) /3, ore = 1 — ff/2. This is line FB in Fig. 03 Above this line, the 
pentagonal region of Fig. 1121 degenerates to a quadrilateral, and the quadrilateral region of Fig. 1131 
degenerates to a triangle. In both cases, vertices (b) and (/) are replaced by a new vertex (/i) 
whose coordinates are a = (3 = (2 — e)/3. These situations are shown in Figs. 1141 and IT51 along 
with corresponding tables of vertex coordinates. 

The quadrilateral of Fig. El w iH degenerate into the triangle of Fig. El when vertices (a) and 
(d) merge and are replaced by vertex (e). As before, this degeneracy happens when e = n, and 
this is the boundary line CJ in the illustration of the (e, 7f) plane shown in Fig. |5J The triangle of 
Fig. 1151 in turn, degenerates into the empty set when (2 — e)/3 = 7f — £, or e = 3ff/2 — 1. Referring 
to Fig. this is the boundary line BC; thus, we see that the quadrilaterals of Fig. El are obtained 
when (e, vf ) £ ACFJ, and the triangles of Fig. El are obtained when (e, 7f) G ABC J. 

Finally, there is one other way that the quadrilaterals of Fig. El can degenerate. Vertices (d) 
and (g) will coincide if (2 — e)/3 = (e — ff)/3, or e = 1 + ff/2. For values of e greater than this, 
vertices (d) and (g) are replaced by vertex (i) with coordinates given by a = —(3 = (2 — e)/3. The 
resulting isoceles triangular region is shown in Fig. 1161 along with a corresponding table of vertex 
coordinates. This triangluar region degenerates to the empty set only when e = 2. Referring to 
Fig. El we see that the triangles of Fig. El are obtained when (e, vf) G ACFL; in fact, since these 
regions are symmetric in (3, they are also symmetric in tt, so their description is the same for all 
(e,7f) G ACDF. 
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Vertex 
(a) 
(b) 
if) 
(<?) 
id) 
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£ + 7f 

3 

2-g 
3 

2-e 
3 

£ — 7f 











£+7T 

3 

37f+4g-2 
6 

37f-4e+2 
6 

— £+7T 



BetaBar 



Figure 12: Bounds on a and /? for 7r = 1/5 and e = 2/3, and the coordinates of the vertices given 
as general functions of e and ff. 
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£ + 7f 
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£+7f 
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2-£ 

3 


3tt+4£-2 
6 


(a) 


2-£ 

3 


3tt-4£+2 
6 


(e) 


7T — 8 


7T — £ 



Figure 13: Bounds on a and (3 for ir = 2/3 and e = 11/20, and the coordinates of the vertices 
given as general functions of e and tt. 
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2-e 
3 


2-e 
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(9) 


2-e 
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37f-4e+2 
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(d) 


e— 7f 
3 


— e+7f 
3 



BetaBar 



Figure 14: Bounds on a and /3 for 7r = 2/3 and e = 1, and the coordinates of the vertices given 
as general functions of e and 7f. 
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2-e 
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2-e 
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37f-4e+2 
6 


(<0 


7f — e 


7T — e 



0.275 0.3 0.325 0.35 0.375 0.4 

BetaBar 



Figure 15: Bounds on a and (3 for -/r = 1 and e = 3/4, and the coordinates of the vertices given 
as general functions of e and 7f. 
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B Basis Bras and Kets for FCHC Model 

In this appendix, we present one possible choice for the bras and kets of the FCHC model for 
three-dimensional fluid dynamics. The basis bras are shown in Table 01 and the basis kets are 
shown in Table |SJ 
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Table 4: Bras for the three-dimensional projection of the FCHC model. The first four rows are the 
hydrodynamic bras, and the last 14 present one choice for the kinetic bras. 
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Table 5: Kets for the three-dimensional projection of the FCHC model, corresponding to the choice 
of bras in Table 0J The first four columns are the hydrodynamic kets, and the last 14 are the kinetic 
kets. 
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